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CONTRACTOR REPORT 


AN IMPROVEMENT IN THE NUMERICAL INTEGRATION PROCEDURE 
USED IN THE NASA MARSHALL ENGINEERING THERMOSPHERE MODEL 


I . INTRODUCTION 


A. A PROBLEM WITH THE MET MODEL 


The models which have been, and shill ace, used to describe 
the properties of the neutral atmosphere between 90 and 
2500km altitude at NASA/Marshall Space Flight Center have 
all been based upon Jacchia's empirical models [1,2]. The 
former model was termed the MSFC/J70 model [3], and in the 
Earth Science and Applications Division of the Structures 
and Dynamics Laboratory the computer program used to output 
data from this model was the J70MM. Recently the computer 
code of the J70MM has been extensively modified, as 
described in [ 4 ] , and the resultant program has been termed 
the NASA Marshall Engineering Thermosphere Model or the MET 
model . 

An intermittent problem in the density output of the 
MSFC/J70 model was found to lie in the integration routine 
which the model employed [4], and this problem was ap- 
parently corrected in the MET model. The MSFC/J70 model and 
the MET model both use Simpson's Rule to numerically in- 
tegrate the barometric and diffusion equations. Although 
the implementation of the integration scheme employed in the 
models is numerically fast, it has not always been found to 
be reliable in the MSFC/J70 model. This unreliability was 
discussed in detail in [4], where it was shown that the con- 
vergence criterion employed was not always stringent enough. 
In an effort to improve reliability, the MET model contained 
a modification in the integration scheme which ensured that 
above a certain altitude a predetermined minimum number of 
iterations were performed before the convergence criterion 
was employed. This modification did improve the reliability 
of the integration method, but further examination has now 
revealed it to be not as reliable as originally thought. 
Therefore a new integration method was sought with the re- 
quirements that it be at least as numerically fast and ac- 
curate as the Simpson's Rule employed in the MET model, but 
that it should be totally reliable. 


One method examined, which is known as Gaussian Quadrature, 
was found to be more than adequate for this purpose. The 
rest of this report is devoted to describing the implementa- 
tion of the Gaussian Quadrature in the MET model as a re- 
placement for the Simpson's integration. Comparisons be- 
tween the two different methods will be made, differences 
will be discussed and finally recommendations will be made. 


B. THE INTEGRAND 


The problem of deciding which numerical method is best to 
use will often depend on the behavior of the integrand one 
is wishing to integrate. Therefore, before discussing the 
integration procedure to be adopted, a short description of 
the behavior of the integrand is in order. 

Between 90 and 105km altitude the density is computed by in- 
tegrating the barometric equation while above 105km altitude 
it is computed by integrating the diffusion equation, as 
described in Jacchia [1]. The barometric equation, valid 
between 90 and 105km altitude is: 


_ z _ 

Mz) - M90) . exp dz} (1) 


while the diffusion equation, valid between 105 and 2500km 
altitude is; 


nj^(z) - nj^(105) 



( 2 ) 


^n Equation (l),p is the mass density, T is the temperature, 
M is the mean molecular weight, g is the gravitational ac- 
celeration and k is the universal gas constant. In Equa- 
tion (2), and M4 are respectively the number density and 
molecular weight of each individual atmospheric specie (N 2 , 
© 2 , O, A, He and H) . Note that in the case of hydrogen 
Equation (2) is applied only above 500km altitude; below 
this altitude njj = 10® m“^ (ie. hydrogen density is 
constant) . 
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One can see that the integrand in Equation (1) is gM/T while 
in Equation (2) it is g/T. The empirical equations for g, T 
and for M below 1051cm altitude, can be found in Jacchia [1]. 
The only adjustable parameters in these equations are the 
exospheric temperature, T^^, and the altitude, z, and thus 
the integrand is only a function of and z. In order to 
examine the behavior of the integrand "two values of were 
chosen which represent, though very approximately, extremes. 
The chosen values were 500K and 2500K. 

The value of gM/T is displayed in figure 1(a) for altitudes 
between 90 and 105 km. The greatest variation of gM/T over 
this altitude range occurs for T„ equal to 2500K, but one 
can see that the variation is quAe small, being typically 
less than 25%. Because the altitude variation of is 
both small and continuous, the behavior of the integrand 
should place no restrictions on the integration method 
employed to integrate the barometric equation. 

The value of g/T is displayed in Figure 1(b) for altitudes 
between 105 and 2500km. Although g/T is a continuous func- 
tion of altitude, its rapid variation at low altitudes may 
tend to make its numerical integration either inefficient or 
inaccurate. If the entire altitude region is separated into 
a number of sub-regions, in such a way that g/T is slowly 
varying in each of these sub-regions, difficulties with the 
numerical integration of g/T will be circumvented. An ex- 
panded representation of this region is shown in Figures 
2(a) and 2(b). By a judicious choice of "sub-regioning, ” 
the altitude variation of g/T will place no restrictions on 
the method employed to integrate the diffusion equation. In 
what follows this is precisely the adopted procedure, but it 
is noted in passing that this procedure has not been adopted 
in the MET model nor the MSFC/J70 model and probably ac- 
counts for some of the problems that these two models have 
experienced when integrating the diffusion equation. 


II. THE GAUSSIAN QUADRATURE 


A. MATHEMATICAL DESCRIPTION 


A concise and useful introduction to Gaussian Quadrature is 
given in [5] and unless explicitly required it is not 
repeated here. 
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ALTITUDE (KM) 



Figure 1(a) The integrand gM/T 
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Gaussian Quadrature can only be applied if integration is 
performed over the interval from -1 to +1, so that generally 
a change of variable is required if this condition is to be 
met. In the problem at hand integration will be between two 
altitudes, say and Z 2 (with Zj<Z 2 ) , so that the following 
transformation will be required: 


Z 


(Za-Zj) 


(x+1) 

2 


+ 


(3) 


Thus, when x=-l, z=z^, and when x=+l, z=Z 2 * Also, 


dz 


1 

2 


(Zj-Zi)dx 


(4) 


Thus, the altitude variable has been changed from z to x. 
In Gaussian Quadrature the transformed integral is ap- 
proximated by the sum of the products of the integrand 
(evaluated at a certain number of discrete points) and cer- 
tain coefficients (associated with each of these points) . 
Mathematically this is stated as 


.+1 N 

f(x) dx « X (5) 

-1 


for some arbitrary function f. Thus, using (4) and (5), 


f 1 ^ 

f(z) dz « X (Zj-Zj) Y. C, f(x ) 

J i-1 


( 6 ) 
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where the integrand f(z) is evaluated at certain discrete 
values of x which correspond to certain values of z in Equa- 
tion (3). The abscissas, Xj^, and the coefficients (or 
weight factors) , C^, for up to eight points (N=8) are given 
in Table 1. It is worth noting that an N-point Gaussian in- 
tegration can evaluate exactly the integral of a 2N-1 degree 
polynomial . 


Table 1. GAUSS QUADRATURE COEFFICIENTS AND ABSCISSAS 


Abscissas = 

= ± : 


Weight 

Factors = Cj^ 

±Xi 


Ci 

±Xi 


Ci 

.57735027 

N*=2 

1.00000000 

.23861919 

.66120939 

.93246951 

N=6 

0.46791393 

0.36076157 

0.17132449 

.00000000 

.77459667 

N=3 

0.88888889 

0.55555556 

.00000000 

.40584515 

.74153119 

.94910791 

N=7 

0.41795918 

0.38183005 

0.27970539 

0.12948497 

.33998104 

.86113631 

N=4 

0.65214515 

0.34785485 

.18343464 

.52553241 

.79666648 

.96028986 

N=8 

0.36268378 

0.31370665 

0.22238103 

0.10122854 

.00000000 

.53846931 

.90617985 

N*5 

0.56888889 

0.47862867 

0.23692689 





B. IMPLEMENTATION IN THE MET MODEL 


The integral of gM/T between 90 and 105km was evaluated very 
accurately using only a 4-point Gaussian integration scheme. 
This is not too surprising if one takes a quick look at 
Figure 1(a), and one also remembers that this scheme will 
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evaluate exactly the integral of a _seventh-degree polyno- 
mial. The mean molecular weight, M, is represented by a 
seventh-degree polynomial in altitude while in this small 
altitude range the gravitational acceleration g, a quad- 
ratic in altitude, is very slowly varying. Temperature T, 
is represented by a fourth-degree polynomial in altitude in 
this altitude range. 

As discussed earlier, the large variations in g/T over the 
altitude range of 105 to 2500km will present a problem for 
most integration methods unless the whole region is 
separated into a number of sub-regions. It was not surpris- 
ing to find, therefore, that even an eight-point Gaussian 
integration (which should exactly evaluate the integral of a 
15-degree polynomial) could not accurately represent this 
integral over the full altitude range, and thus a number of 
sub-regions were constructed. 

Within each sub-region the integral of g/T was evaluated by 
an N-point Gaussian integration scheme, where the number of 
points used in the integration, N, was determined by ex- 
perimentation for each particular sub-region. After much 
experimentation the whole region extending from 105 to 
2500km was separated into seven sub-regions. The altitude 
extent of each of these regions as well as the number of 
points, N, used in the Gaussian integration for each region 
is given in Table 2. For completeness, the region below 
105km is also included. 


Table 2. DEFINITION OF SUB-REGIONS 


Altitude Range 
(km) 


90- 105- 125- 160- 200- 300- 500- 1500- 
150 125 160 200 300 500 1500 2500 


No. of points, 
N, in Gaussian 
integration 


4566666 6 


The incorporation of the Gaussian Quadrature into the MET 
model was affected primarily by replacing subroutine IN- 
TEGRATE with a new subroutine named GAUSS. This new sub- 
routine handles the setting up of the sub-intervals (ie. the 
altitude boundaries and the number of points required for 
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the Gaussian Quadrature for each sub-interval ) , and then 
performs the Gaussian Quadrature in only seventeen lines of 
executable code! Subroutine GAUSS can be found in the Ap- 
pendix along with the associated FORTRAN coding of the MET 
model . 


Ill . RESULTS 


In this section the integral and total mass densities are 
evaluated by the standard integration method employed in the 
MET model (using Simpson's Rule) and by the new method 
(Gaussian Quadrature) so that the two different integration 
methods can be compared. For these two methods, the ac- 
curacy and reliability of the results and the computational 
speeds will be compared. To ensure that the comparison is a 
reliable one, it would be advantageous to know the exact 
value of the desired quantity (the value of an integral or 
the total mass density) . Since this is not known, a 
"reference” version of the MET model was written in double- 
precision. This reference model employed Simpson's Rule in 
the same way that the MET model does, but, with three excep- 
tions ._ Firstly, the convergence parameter (c) is set equal 
to 10 in the MET model, but in the double-precision 
reference model, it is set to 10”" (it could have been set 
even smaller than this, but the model would then have con- 
sumed a prohibitive amount of CPU time) . Secondly, in the 
reference model, once convergence has been achieved, one 
more iteration is performed and the convergence condition is 
tested a second time. Two consecutive positive convergence 
tests constitute convergence in this model. Lastly, whereas 
the maximum number of sub-intervals used in the Simpson's 
integration in the MET model is 2^®, the maximum number that 
can be used in the reference model is 2^^. These changes 
ensure that the results obtained from the reference model 
will be both accurate and reliable, but it should be remem- 
bered that this is achieved at the expense of CPU time. 

To present results for this report, two different sets of 
input parameters were used in the models. The first set of 
input parameters, which lead to a cold atmosphere with an 
exospheric temperature of 597.361K, is as follows: 
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Date & Time: 

^10.7 ^ ^10.7* 
Latitude : 


1987, June 21, 0400 hr UT 
70 Ap index: 0 

0 Longitude : 0 


The second set of input parameters lead to a hot atmosphere 
with an expspheric temperature of 2554. 217K and is given as: 


Date & Time: 

^10.7* 

Latitude: 


1987, October 27, 1400 hr UT 
400 ^10. 7* ^p iudex: 
0 Longitude: 0 


400 


Although the exospheric temperatures for the cold and the 
hot atmospheres are not the same as those which were used to 
represent the extreme conditions shown in Figures 1 and 2, 
the general behavior of the integrands will not be sig- 
nificantly different so that the previous discussion of the 
integrands will still apply. A good description of these 
input parameters can be found in [3], and will not be 
repeated here. 

All of the results presented here (excepting those in sub- 
section E) are represented graphically in figures 3 through 
to 10. In order to maximize the sensitivity of the com- 
parisons, a large number of data points were used in all of 
the figures. In the 90-105km section nearly 100 data points 
were generated while in the 105-2500km section some 800 data 
points were generated. 


A. COMPARISON OF THE INTEGRALS: THE 90 TO 105KM SECTION 


In this altitude range the integral of gM/T is evaluated by 
the three methods already discussed. The value of this in- 
tegral calculated using the standard Simpson's Rule as used 
in the MET model is designated RS, and that calculated using 
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the reference model (double-precision, high accuracy 
Simpson's Rule) is designated RD, while that calculated 
using Gaussian Quadrature is designated RG. This terminol- 
ogy will be used throughout this report. 


In order to compare RS and RG, the percentage deviations of 
each of these from RD are calculated as: 


% deviation = (R/RD - 1) x 100% 


(7) 


where R is equal to either RS or RG. Later in this report, 
density deviations will be calculated in a similar fashion. 
Assuming the high-accuracy reference model gives correct 
results means that these percentage deviations can then be 
considered as errors in the other two integration methods, 
and henceforth they will be referred to as such. 

The error in RS as a function of altitude in the cold atmos- 
phere is shown in Figure 3(a). At low altitudes, i:he errors 
in RS are small, typically less than about 3x10“^% in mag- 
nitude, but they increase steadily with increasing altitude 
and approach a value close to 10“^%, just below 102km al- 
titude. Just after the maximum error is reached, the errors 
decrease rapidly with an increase of altitude, and the 
process appears to repeat itself again. The reason for this 
behavior is quite simple. The number of sub-intervals used 
in the Simpson's Rule, which calculates RS, is the same for 
all integration intervals extending from 90km to anywhere at 
or below the altitude where the maximum error is achieved. 

Above that altitude the number of sub-intervals is doubled. 
The error is thus drastically reduced, but as one progresses 
to greater altitudes, the error begins to increase again. 
This behavior will be observed repeatedly in all of the 
results employing Simpson's Rule which follow, both in RS 
and in density. 

The error in RG as a function of altitude in the cold atmos- 
phere is shown in Figure 3 (b) . The errors in RG are small 
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at all altitudes, never exceeding about 3x10“®% in mag- 
nitude. Most of this error is probably associated with 
single-precision round-off. because on the VAX 11/780 a 
precision of one part in 10' is typically to be expected in 
the representation of a single-precision real number. A 
comparison of figures 3(a) and 3(b) reveals that at very low 
altitudes (between 90 and 95km altitude) the errors in using 
Simpson's Rule are very similar to those arising from the 
use of the Gaussian Quadrature. However, at greater al- 
titudes, the errors involved in the use of Simpson's Rule 
increase remarkably, while those associated with the use of 
the Gaussian Quadrature remain essentially small and 
bounded . 

The results corresponding to the hot atmosphere are shown in 
Figures 4(a) and 4(b). There are similarities between these 
results and those shown in Figures 3(a) and 3(b), but once 
again the overall trend is that the use of the Gaussian 
Quadrature gives much smaller errors than through the use of 
Simpson's Rule. Typically more than an order of magnitude 
smaller. 

Extensive timing tests revealed that the Gaussian integra- 
tion was always faster than the Simpson's integration. The 
integrand was always evaluated four times in the Gaussian 
integration. In the Simpson's' integration the integrand 
was evaluated four times in the altitudes below where the 
maximum errors occurred (-100km) , and six times for al- 
titudes above where the maximum errors occurred. At 105km 
altitude the integration took twice the CPU time using 
Simpson's Rule, but the errors in RS were about an order of 
magnitude larger than those in RG. Hence, in the 90-105km 
altitude range the Gaussian integration can be twice as fast 
and ten times more accurate than the Simpson's integration. 


B. COMPARISON OF THE INTEGRALS: THE 105 TO 2500KM SECTION 


In the 105 to 2500km altitude range the integral of g/T is 
evaluated by the three methods. The error in RS as a func- 
tion of altitude in the cold atmosphere is shown in Figure 
5(a). One notices that at 416km and around 1350km altitude 
the errors in RS are extremely large, and are much larger 
than expected from the size of the convergence parameter 
(10“^) . These errors are not associated with round-off, but 
are due to the unreliability of the Simpson's' integration 
method as implemented in the current MSFC/J70 and NASA MET 
models. 
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Figure 5(a) 
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Figure 5. Error in the integral of g/T evaluated using (a) and (b) Simpson's Rule and 
(c) Gaussian Quadrature for the cold atmosphere. See text for details. 
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The unreliability of this method, which is due to an oc- 
casional and unpredictable false convergence condition in 
the iteration scheme, is discussed in [4], At other al- 
titudes, however, the Simpson's integration method appears 
to be reliable. In order to examine the errors associated 
with the reliable results, the erroneous results were 
changed to the same values as those of the neighboring al- 
titudes. These errors, which were undiscernable in figure 
5(a), are shown in Figure 5(b). The average errors appear 
to have a magnitude of 1 to 2 x 10“^%, but maximum errors 
can be as large as 8 x 10“^% in magnitude. Once the number 
of sub-intervals used is doubled, tlje errors reduce drasti- 
cally to magnitudes of about 3 x 10“’®%. 


The error in RG, as a function of altitude, in the cold at- 
mosphere is shown in Figure 5(c). One notices that the er- 
rors associated with the Gaussian integration remain small 
at all altitudes, being on average less than 10”®% in mag- 
nitude. The ^argest errors which occur are a little larger 
than 2 X 10“®%, which is probably due mostly to round-off 
error. The Gaussian integration method appears to be to- 
tally reliable, showing no evidence of the type of errors 
which can occur with integration by Simpson's Rule (see 
Figure 5(a)). When the integration using Simpson's Rule is 
successful, the errors associated with it are on average an 
order of magnitude larger than those associated with the 
Gaussian integration. 


The equivalent set of results for the hot atmosphere are 
shown in Figures 6(a), (b) and (c) . Under these conditions 
the Simpson's integration appears to be more unreliable than 
it was for the cold atmosphere (Figure 6(a)); and once again 
integration by Simpson's Rule introduced errors which are on 
average an order of magnitude larger than these associated 
with the Gaussian integration. 


Extensive tests which measured the amount of CPU time that 
was required to perform the integrations revealed that the 
Gaussian integration was always faster than the Simpson's 
integration. It was faster by a factor of approximately 2 
or 3 at 200km altitude (2 in the hot atmosphere and 3 in the 
cold atmosphere), approximately 4.5 at 500km altitude, 7.4 
at 1000km altitude increasing to about 12.5 at 2500km al- 
titude. Thus, integration using Gaussian Quadrature will 
save a considerable amount of CPU time over the use of 
Simpson's Rule. 
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Figure 6(a) 


Figure 6(b) 


Figure 6(c) 
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Figure 6. Error in the integral of g/T evaluated using (a) and (b) Simpson’s Rule and 
(c) Gaussian Quadrature for the hot atmosphere. See text for details. 






In the following two sub-sections, the densities calculated 
for the two altitude regions using Gaussian Quadrature and 
using Simpson's Rule are compared in a similar way that the 
comparisons were performed in the previous section. 


C. COMPARISON OF THE DENSITIES: THE 90 TO 105KM SECTION 


The errors In the density values calculated using Simpson's 
Rule for the cold atmosphere between 90 and lOSJcm altitude 
are shown in Figure 7(a), while those associated with the 
Gaussian Quadrature are shown in figure 7 (b) . The magnitude 
of the maximum error associated with the use of Simpson's 
Rule is almost 2 x 10“^%. Over almost half of the altitude 
range the magnitude of the errors associated ^th the use of 
Simpson's Rule is greater than about 2 x 10“^%, whereas in 
the case of the Gaussian Quadrature these errors are ’’zero” 
(within the accuracy of single-precision arithmetic) over 
most (all but 2 km) of the altitude range. For the hot at- 
mosphere, very similar results are obtained (see Figures 
8(a) and 8(b)). Thus, independent of atmospheric conditions 
(hot or cold) , the densities are calculated significantly 
more accurately using Gaussian Quadrature to integrate the 
barometric equation than by using Simpson's Rule in the 90- 
105km altitude region. 


D. COMPARISON OF THE DENSITIES: THE 105-2500KM SECTION 


The errors in the density values calculated using Simpson's 
Rule for the cold atmosphere between 105 and 2500km altitude 
are shown in Figure 9(a). The large errors at 416km and 
1350km altitude are associated with the unreliability of the 
Simpson's integration method which gave large errors in the 
calculations of the integral of g/T, as previously shown in 
Figure 5(a). Once again, in order to examine the errors as- 
sociated with reliable integration using simpson's Rule the 
large errors were removed from Figure 9(a), and the remain- 
ing errors in the density values were re-plotted in Figure 
9(b). The magnitudes of the errors are usually smaller than 
3 X 10“^ % except just above 400km altitude where they rise 
to almost 5.5 X 10“^%. 

The corresponding errors in the density values calculated 
using Gaussian Quadrature are shown in figure 9(c). These 
errors are negative at all altitudes, meaning that the cal- 
culated densities are smaller than they should be. The 
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Figure 7. Error in the density evaluated using (a) Simpson’s Rule and (b) Gaussian Quadrature 
for the cold atmosphere below lOSkm altitude. 
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Figure 8. 
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Error in the density evaluated using (a) Simpsons’s Rule and (b) Gaussian Quadrature 
for the hot atmosphere below lOSkm altitude. 
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Figure 9. Error in the density evaluated using (a) and (b) Simpson’s Rule (c) Gaussian 
Quadrature for the cold atmoshere above 105km altitude. See text for details. 






magnitudes of the errors also generally increase with in- 
creasing altitude, being less than about 2.5 x at low 
altitudes and rising to values of about 3 x 10 at the 
highest altitudes. A comparison of Figures 9(c) and 9(b) 
reveals that at the highest altitudes (say, greater than 
1000km) the errors arising from the two integration methods 
(Simpson's Rule and Gaussian Quadrature) are much the same, 
providing that one overlooks the errors of about 3 x 10~^% 
centered around 1350km altitude associated with the un- 
reliability of the Simpson's integration method employed in 
the MET model. At low altitudes, however, the densities 
calculated using the Gaussian Quadrature will always be ac- 
curate, usually much more accurate than those calculated 
using Simpson's Rule. 

The corresponding set of results for the hot atmosphere are 
shown in Figures 10(a), (b) and (c) . In the hot atmosphere 
the Simpson's integration method appears to be considerably 
more unreliable than it is in the cold atmosphere (compare 
Figures 9(a) and 10(a)). A comparison of figures 10(c) and 
9(c) shows that the Gaussian Quadrature gives much more ac- 
curate density values in the hot atmosphere than in the cold 
atmosphere. In the hot atmosphere the magnitudes of these 
errors are usually less than 4 x 10~^%, so that the Gaussian 
Quadrature calculates densities significantly more ac- 
curately than does the Simpson's integration method at all 
altitudes in the hot atmosphere. 


E. A NUMERICAL EXAMPLE 


Here, the density value output from the standard MET model 
using Simpson's Rule, modified MET model using Gaus- 

sian Quadrature, Pq , and the high-precision reference model 
(double precision, smaller tolerance, etc.), , are given. 
These results are for the hot atmosphere as defined by the 
input parameters given on page 10, and for an altitude of 
400km: 

Pg - 3.3844014 X lO'^^kgm"^ 

Pq “ 3.3844164 x 10 ^^kgm ^ 
p^ - 3.3844164 x lO'^^kgm'^ 
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Figure 10(c) 



Figure 10. Error in the diensity evaluated using (a) and (b) Simpson’s Rule and (c) Gaussian 
Quadrature for the hot atmosphere above 105km altitude. See text for details. 






Thus, the error in the density value calculated with the MET 
model is about -4.43 x 10”^% while the corresponding error 
with the modified MET model which uses Gaussian Quadrature 
is zero. 


IV. DISCUSSION 

The evaluation of the integral of gM/T between 90 and 105km 
altitude and g/T between 105 and 2500km altitude has been 
performed by two different methods. The first of these 
methods employs a modified form of Eimpson's Rule, basically 
as used in the MSFC/J70 and NASA MET models. The second 
method is Gaussian Quadrature. Extensive comparisons of the 
results of using these two methods have revealed that in- 
tegration by Gaussian Quadrature is significantly numeri- 
cally faster than by Simpson's Rule and that the integration 
by Gaussian Quadrature is overall significantly more ac- 
curate than by Simpson's Rule. It was also found that at 
certain altitudes the form of Simpson's Rule used in the MET 
model becomes unreliable, where errors in the results become 
significantly larger than expected from the convergence 
criterion used in the iterative scheme of the method. This 
unreliability occurred more often in the hot atmosphere. No 
such unreliability was evident in any of the results which 
involved Gaussian Quadrature. 

Extensive comparisons of the densities calculated using the 
two methods revealed, not surprisingly, that those resulting 
from the use of Gaussian Quadrature were overall sig- 
nificantly more accurate than those resulting from the use 
of Simpson's Rule. Also, reliability of the calculated den- 
sities was a problem with Simpson's Rule, but not with the 
Gaussian Quadrature. Any program which is written to calcu- 
late densities will run faster if integration is done using 
Gaussian Quadrature, however, just how much faster it will 
run, than an equivalent program which integrates using 
Simpson's Rule, will depend on the nature of the particular 
program involved. When Gaussian Quadrature was used, the 
program which was used to generate the data files, which are 
plotted in Figures 7(a) through to 10(c), ran only very 
slightly faster at low altitudes and ran five times faster 
above 2000km altitude. 

To all intents and purposes the results obtained by using 
the Gaussian Quadrature are as accurate and reliable as 
those obtained from the "reference” model. One must recall 
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that this "reference" model employed the same form of 
Simpson's Rule that is used in the MET model, but it is 
written in double-precision. It has a convergence parameter 
that is one-thousand times smaller than that used in the MET 
model, and requires completion of two consecutive positive 
convergence tests before convergence is accepted. It is no 
wonder that integration using this reliable version of 
Simpson's Rule took anywhere from about 30 to 400 times 
longer than by using Gaussian Quadrature. 

One should also note that the choice of sub-interval 
parameters which are given in Table 2 could possibly be im- 
proved upon or modified to suit ones individual needs. For 
general use, however, the values given in Table 2 are recom- 
mended because they do not compromise on accuracy, and the 
Gaussian Quadrature using these values still requires less 
computational time than the standard method (Simpson's Rule) 
which is currently used in the MSFC/J70 and MET models. 

Finally, it should be noted that if integration is performed 
between two different altitudes (say Zi, and Z 2 ) and the 
fractional errors associated with thi^ integration are 
greater at one altitude (z,) than the other (Z 2 ) , it does 
not necessarily follow that when the densities are calcu- 
lated the fractional error in the density at altitude z^ 
will be greater than that at Z 2 * This is because the den- 
sity calculation involves integration errors being multi- 
plied by numbers which are altitude dependent. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Although the reliability of the integration scheme used in 
the MSFC/J70 model was improved in the MET model, the 
results presented in this report have shown that it is still 
not as reliable as desired. An alternative integration 
scheme, based on Gaussian Quadrature, has been substituted 
for the standard integration scheme used in the MET model 
and has been found to be totally reliable. At low al- 
titudes, it generally gives significantly more accurate den- 
sity results than does the standard MET model, and uses less 
computer time to do it. At high altitudes, it generally 
gives results which have a similar accuracy to those of the 
standard MET model, but it uses only a fraction of the com- 
puter time to do it. 

Due to the reliability, the accuracy, and the speed of this 
new integration scheme, which uses Gaussian Quadrature, it 
is recommended that this new scheme replace the older, more 
unreliable scheme, based upon Simpson's Rule, as presently 
used in the MET model. If this is done, it is emphasized 
that the integration parameters given in Table 2 should be 
adhered to. This being the case, it is further recommended 
that the program listed in the Appendix, which is also named 
the NASA Marshall Engineering Thermosphere (MET) Model, and 
which contains the essentials of this new integration scheme 
in subroutine GAUSS, should replace the current MET program. 
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APPENDIX 


c* 

C* The Marshall Space Flight Center 

C* Marshall Engineering Thermosphere Model 

C* 

c* 

C* written by 

Ce 

C* Mike Hickey 

C* Universities Space Reseorch Association 

C* NASA / MSFC . ED44 

Ce TeK (205) 544-5692 

C* 

C* This program is a driving program for the following subroutines 

C* 

C* ATMOSPHERES 

C* SOLSET 

C* TIME 

C* J70 

C* 

C* The atmospheric model is a modified Jacchio 1970 model and is given in 

C* the subroutine J70. All of the other subroutines were designed to 

C* allow flexible use of this model so that various input porometers could 
C* be varied within a driving program with very little software development. 

C* Thus, for example, driving routines can be written quite easily to 

C* facilitote the plotting of output os line or contour plots. Control is 

C* achieved by setting the values of four switches in the driving program, 

C* os described in subroutine ATMOSPHERES. 

C* 

C**«*****«***«*********«««*«*«**««*e«»«»**********«*«*«***e*****e*e*********** 


REAL*4 INDATA (12) , OUTDATA (12) , AUXDATA (5) 

CHARACTERel SWITCH (4) 

CALL LIB$INIT_TIMER 

C Set all switches to *Y* so thot only one particular calculotion is performed 

SWITCH (1) - 'Y* 

SWITCH (2) - *Y' 

SWITCH (3) - •Y* 

SWITCH (4) - *Y* 

CALL ATMOSPHERES ( INDATA, OUTDATA. AUXDATA. SWITCH ) 

C Now type output date 


Typ* 


All output In MKS units* 




Type 

•.* 

• 





Type 


Exospheric temperature 

- 

• 

f 

OUTDATA 

(1).* K* 

Type 


Temperature 

■1 

» 

t 

OUTDATA 

(2).' K* 

Type 


N2 number density 

m 

• 

f 

OUTDATA 

(3).' /m3* 

Type 


02 number density 


• 

• 

OUTDATA 

(4).* /i»3* 

Type 

*.• 

0 number density 

m 

• 

• 

OUTDATA 

(5),* /m3‘ 

Type 


A number density 

m 

• 

• 

OUTDATA 

(6).' /n.3' 

Type 

•.* 

He number density 

m 

• 

• 

OUTDATA 

(7).' /m3* 

Type 


H number density 

H 

f 

t 

OUTDATA 

(8).* /m3* 

Type 


Average mol ecu lor wt. 

■i 

• 

• 

OUTDATA 

(9) 

Type 


Total mass dens i ty 

m 

* 

• 

OUTDATA 

(ie),’ kg/iii3* 

Type 


Log10 mass density 


• 

• 

OUTDATA 

(11) 

Type 


Total pressure 


• 

• 

OUTDATA 

(12),* Pa’ 

Type 


Local grav. acceln. 


t 

> 

AUXDATA 

(1).' m. sec-2’ 

Type 


Ratio specific heats 

m 

• 

• 

AUXDATA 

(2) 

Type 


Pressure scale-height 

B 

• 

• 

AUXDATA 

(3).’ m’ 

Type 


Specific heat cons, p 

B 

• 

• 

AUXDATA 

(4) , ’ib2.8SC-2.K— 1 

Type 


Specific heat cons, v 

B 

• 

» 

AUXDATA 

(5),’m2.ssc-2.K-1 


Type • 



call LiB$skmjnum 


STOP 

END 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


SUBROUTINE ATMOSPHERES ( INDATA, OUTDATA, AUXDATA. SWITCH ) 


DESCRIPTION:- 


Calculate atmospheric data in single precision using subroutine 070 
ond J70SUP. 


SUBROUTINES: 


TIME. SOLSET, CMC, J70 and J70SUP 
INPUT:- 


all single precision, either through 
subroutines or from main driver prog. 


INDATA (1) — altitude - Z 

(2) — latitude - XLAT 

(3) — longitude * XLNG 

U) — y*or (yy) - iyr 

(5) — month (mm) * MN 

(e) — day (dd) - IDA 

(7) — hour (hh) - IHR 

(8) — mins (mm) ■ MIN 

(9) — geomognetic index ■ IGEO^IND 

(10) — solar radio noise flux» F10 

(11) — 162-day average F10 « F10B 


(12) — geomagnetic activity index « GI-^ 


OUTPUT :- 


NOTE : All output in MKS units 


al I single precision 


OUTDATA (1) — exospheric temperature (K) 

(2) — temperature at altitude Z 

(3) — N2 number density (per meter-cubed) 


(4) — 02 number density ( .• ) 

(5) — 0 number density ( .. ) 

(6) — A number density ( .. ) 

(7) — He number density ( .. ) 

(8) — H number density ( .. ) 

(9) — average molecular weight 

(10) — total density 

(11) — Iog10 ( total density ) 

(12) — total pressure ( Po ) 


AUXDATA (1) — gravitational acceleration ( m/s-s ) 

(2) — ratio of specific heats 

(3) — pressure scale-height ( m ) 

(4) — specific heat at constant pressure 

(5) — specific heat at constont volume 


COMylENTS:- 


SWITCH(I) — if Y(es), 
SWITCH(2) — if Y(es), 
SWITCH(3) — if Y(es), 
SWITCH(4) — if Y(es), 


date and time are input from terminal through 
subroutine TIME once only 

solor/mognet ic activity are input from terminal 
through subroutine SOLSET once only 
only ONE altitude value is input from terminal 
through main calling program 

only ONE latitude AND longitude ore input from 
terminal through main calling program 


30 



oo oo oo oo oo 




C* ATMOSPHERES written by Mike Hickey ( USRA, NASA/ED44 ) 
C* Tel: (205) 544-5692 

C* _____ Jonuory-Aprl 1 1987 




EXTERNAL TIME 
DIMENSION AUXDATA (5) 

INTEGER HR 

REAL«4 LAT. LON, INDATA (12). OUTDATA (12) 
CHARACTERel SWITCH (4) 

PARAMETER PI - 3.14159265 


Thie next eectlon Is only executed on the first coll to ATMOSPHERES 
DO WHILE ( CALL. EQ. 0.0 ) 


SECTION A: 


IF ( SWITCH(I). EQ. 'Y* ) THEN 

CALL TIME ( lYR, MON. IDA. HR, MIN. SWITCH(I) ) 
INDATA (4) - FLOATJ UYR) 

INDATA (5) - FLOATJ (MON) 

INDATA (6) - FLOATJ (IDA) 

INDATA (7) - FLOATJ (HR) 

INDATA (8) - FLOATJ (MIN) 

END IF 


SB^TION B:- 


IF ( SWITCH(2). ra. *Y' ) THEN 

CALL SOLSET ( IGEO.IND. F10. F10B, GI. SWITCH(2) ) 
INDATA (9) - FLOATJ (IGEO.IND) 

U«>ATA (10) •* F10 
INDATA 01) - H0B 
INDATA (12) - GI 

END IF 


SECTION C:- 


IF ( SWITCH(3). EQ. *Y* ) THEN 

TYPE *.* Input oltltude, las* 
ACCEPT e. INDATA (1) 

Z • INDATA (1) 

END IF 

SECTION D:- 


IF ( SWITCH(4). EQ. *Y* ) THEN 

TYPE •.* Input lotitude and lon 9 ltude, degrees* 

ACCEPT *. ( INDATA(I), !• 2.3 ) 

LAT • INDATA (2) 

LON - INDATA (3) 

RLT “ INDATA (2) e PI / 180. I geogrophle lotitude, rodlons 
END IF 
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CALL - 1 .• 

END DO 

End of first executable section 


C The following depend on the values of the switches 


C*«** 

C* SECTION 1:- 


IF ( SWITCH(I). NE. *Y* ) THEN 

lYR - JINT ( INDATA (4) ) 

MON - JINT ( INDATA (5) ) 

IDA - JINT ( INDATA (6) ) 

HR - JINT ( INDATA (7) ) 

MIN - JINT ( INDATA (8) ) 

CALL TIME ( lYR, MON, IDA. HR, MIN. SWITCH(I) ) 
END IF 


C**«** 

C* SECTION 2:- 


IF ( SWITCH(2). NE. *Y* ) THEN 

IGEO.IND - JINT ( INDATA (9) ) 

Fie - INDATA (ie) 

F1CB - INDATA (11) 

GI - INDATA (12) 

CALL SOLSET ( IGEO.IND, F10, F10B. GI, SWITCH(2) ) 
END IF 


C*«*** 

C* SECTION 3:- 


IF ( SWITCH(3). NE. *Y* ) THEN 
Z - INDATA (1) 

END IF 


C**«** 

C* SECTION 4:- 


IF ( SWITCH(4). NE. ’Y* ) THD4 

UT - INDATA (2) 

LON - INDATA (3) 

RLT > INDATA (2) * PI / 180. I geographic lotitude, radians 
END IF 


C All setting-up complete. 

CALL J70 ( INDATA, OUTDATA ) 

CALL J70SUP ( 2. OUTDATA. AUXDATA ) 


RETURN 

ENTRY ATMOS.ENT ( DIMIY ) 

CALL - DUMMY 

RETURN 

END 
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SUBROUTINE TIME ( lYR, MON. IDA. HR. MIN, SWITCH ) 


C- 


c* 

c* 


This subroutine sets up time of year and day 

c* 

c* 


INPOTS/OUTPUTS: 

c* 

c* 

lYR 

« year ( 2 digits ) 

c* 

MON 

^ month 

c* 

IDA 

— day of month 

c* 

HR 

- hour of day 

c* 

MIN 

■ minutes 

c* 

c* 


Written by Mike Hickey, USRA 

c** 



DIMENSION IDAY ( 12 ) 

INTEGER HR 
CHARACTER* 1 SWITCH 

DATA IDAY / 31. 28. 31. 38. 31. 38 . 31, 31. 38 . 31, 38 , 31 / 
PARAMETER PI - 3.14159265 

O- — 

C If SWITCH - Y(e«) then Input data and time from terminal 

C 


IF ( SWITCH. EQ.’Y*. OR. SWITCH. EQ. ) THEN 

TYPE *. * Input date and time of date? ( yy.mm.dd.hh.mm ) * 
ACCEPT *. lYR, MON. IDA. HR. MIN 

END IF 


IF ( JMOD (IYR.4) .EQ. 8 ) THEN 

IF ( JMOD (IYR.188) .NE. 8 ) IDAY ( 2 ) - 29 

ELSE 

IDAY ( 2 ) - 28 

END IF 

DAYTOT - 8.8 


DO 1 I - 1. 12 

DAYTOT - DAYTOT + FLOATJ ( IDAY ( I ) ) 

1 CONTINUE 


IF ( MON. GT. 1 ) THEN 


KE « MON - 1 
ID i- 8 

DO 2 I « 1 . KE 
ID - ID + IDAY (I) 

2 CONTINUE 


ID -i ID 4> IDA 
DO • IDA 

END IF 


ELSE 


RETURN 

END 
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SUBROUTINE SOLSET ( IGEO^IND. F10, F10B, GI, SWITCH ) 


c* 

C* This subroutine simply colls for o setup of the solor-octivty and auroral 
C* activity indices. 


c* 

c* 


INPUTS/OUTPUTS; 

c* 

c* 

IGEO.IND - 

geomagnetic index 

c* 

F1C 

solar radio noise flux 

c* 

F10B 

162-day average F10 

c* 

GI 

geomagnetic activity Index 

c* 



c* 

c** 

Written by Mike Hickey, USRA 


CHARACTER*! SWITCH 
IGEO^IND - 2 

C 

C If SWITCH - Y(es) then Input geomagnetic indices from terminal 

C, 

IF ( SWITCH. EQ.’V. OR. SWITCH. EQ. ’y* ) THEN 

C TYPE *, • Input geomagnetic Index ( 1-KP, 2-AP ) * 

C ACCEPT *. IGEO^IND 

TYPE *, • Input solar radio noise flux ( F10 - 0-400 ) * 

ACCEPT *. F10 

TYPE *. • Input 162-day average F10 ( F10B - 0-250 ) ' 

ACCEPT *, F10B 

C IF ( IGEO^IND . EQ. 2 ) THEN 

C 

C TYPE *, • Input geomagnetic octivlty index ( GI « 0-400 ) • 

C 

C ELSE 

C 

C TYPE *, • Input geomagnetic activity index ( GI - 0-9 ) * 

C 

C END IF 

TYPE *,• Input AP Index ( AP - 0 - 400 ) • 

ACCEPT *. GI 

END IF 

C ! 


RETURN 

END 



SUBROUTINE J70SUP ( 2, OUTDATA, AUXDATA ) 


C** 

C* 

C* 

C* 

C* 

C* 

c* 

c* 

c* 

€♦ 

C« 


DESCRIPTION:- 


J70SUP calculates auxilliary variables which are output in orray 
AUXDATA, given data input from J70 which are contained in array OUTDATA. 

INPUT DATA:- 


C* 

2 

— 

altitude (km) 


OUTDATA (2) 

c* 

T22 

— 

temperature at altitude z 

- 

c* 


— 

N2 number density 

■1 

.. (3) 

c* 


— 

02 . . 

■ 

.. (4) 

c* 


— 

0 .. 

- 

.. (5) 

c* 


— 

A .. 

- 

.. (6) 

c* 


— 

He . . 


.. (7) 

c* 


— 

H .. 

Wt 

.. (8) 

c* 

EM 

— 

average moteculor weight 

- 

.. (9) 

c* 

DENS 

— 

total density 

- 

.. (ie) 

c* 

P 

— 

total pressure 

- 

.. (12) 


Ce 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

0 * 

C* 


OUTPUT DATA:- 


G — gravitational acceleration » AUXDATA (1) 

GAM — ratio of specific heats ■ AUXDATA (2) 

H — pressure scale-height ■ AUXDATA (3) 

CP — specific heat at constant pressure « AUXDATA (4) 

CV — specific heat at constant volume ■ AUXDATA (5) 


Written by Mike Hickey^ USRA 


REAL*4 OUTDATA (12). AUXDATA (5), H 

G - 9.B0665 / ( ( 1. + 2 / 6.356766E3 )**2 ) 

H - OUTDATA (12) / ( G • OUTDATA (10) ) 

SUM1 « OUTDATA (3) + OUTDATA (4) 

SUM2 - 0.0 

DO 1 I - 5, 8 
SUM2 - SUM2 + OUTDATA (I) 

1 CONTINUE 

gam - ( 1.4 • SUM1 4 1.67 e SUM2 ) / ( SUM1 4 SUM2 ) 

CV-GeH/( (GAM-I.e)e OUTDATA (2) ) 

CP - GAM * CV 

AUXDATA (1) « G 
AUXDATA (2) - GAM 
AUXDATA (3) - H 
AUXDATA (4) - CP 
AUXDATA (5) - CV 

RETURN 

END 
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SUBROUTINE J70 ( INDATA, OUTDATA ) 


c** 



• • 

c** 

J70 developed from J70MM by 

mm 

c** 

Mike P. Hickey 


mm 

c** 

Universities Space Reseorch 

Assoc i at i on 

mm 

C*e 

at 


mm 

C*e 

NASA / MorshafI Space Flight 

Center, ED44, 

mm 

C** 

HuntsviltSe Aloboma, 35812, USA. 

mm 

C** 

Tel. (205) 544-5692 


mm 

C*e 



mm 

C** INPUTS: 

through the subroutine cal 

ling list 

mm 

C** 



mm 

C*» OUTPUTS: 

through the subroutine cal 

ling list 

mm 

c** 



mm 

c** 



mm 

0** 

INPUT DATA: 


mm 

c** 



mm 

c** 

2 — altitude 

- INDATA (1) 

mm 

c** 

XLAT — latitude 

- INDATA (2) 

mm 

c** 

XLNG — longitude 

- INDATA (3) 

mm 


lYR — year (yy) 

- INDATA (4) 

mm 

c** 

MN — month (mm) 

- INDATA (5) 

mm 

c** 

IDA — day (dd) 

- INDATA (6) 

mm 


IHR — hour (hh) 

- INDATA (7) 

mm 

0** 

MIN — mins (mm) 

- INDATA (8) 

mm 

c** 

11 — geomagnetic index 

- INDATA (9) 

mm 

c** 

F10 — solar radio noise flux 

- INDATA (10) 

mm 


F10B — 162-day average F10 

- INDATA (11) 

mm 

c** 

GI — geomagnetic activity index - INDATA (12) 

mm 

C*e 



mm 

C** 



mm 


OUTPUT DATA: 


mm 

C^e 



mm 


T — exospheric temperature 

- OUTDATA (1) 

mm 

C** 

TZZ — temperature at altitude Z 

- OUTDATA (2) 

mm 

C** 

A(1) — N2 number density 

- OUTDATA (3) 

mm 


A(2) — 02 number density 

- OUTDATA (4) 

mm 


A(3) — 0 number density 

- OUTDATA (5) 

mm 


A(4) — A number density 

- OUTDATA (6) 

mm 


A(5) — He number density 

- OUTDATA (7) 

mm 

0** 

A(6) — H number density 

- OUTDATA (8) 

mm 

c** 

DA — averoge molecular weight 

- OUTDATA (9) 

mm 


DENS — total density 

- OUTDATA (10) 

mm 

C** 

DL — Iog10 ( total density ) 

- OUTDATA 01) 

mm 

C** 

P — total pressure 

- OUTDATA (12) 

mm 

C** 



mm 

NB. 

Input through or ray ’INDATA* 


mm 


Output through array ’OUTDATA* 


mm 


DIMENSION A ( 6 } 

REAL*4 INDATA ( 12 ). OUTDATA ( 12 ) 

PARAMETER RGAS - 8.31432E3 I J/kmol-K 

PARAMETER BFH - 440.C 

C Calcultions performed for only one latitude , one longitude 
C and one oltitude 


C 


Cmm Set parameters to 
r 

INDATA 

values 

z 

■I 

INDATA 

(1) 



XLAT 

m 

INDATA 

(2) 



XLNG 

m 

INDATA 

(3) 



lYR 

m 

JINT ( 

INDATA 

(4) 

) + 1900 

MN 

m 

JINT ( 

INDATA 

(5) 

) 

IDA 

m 

JINT ( 

INDATA 

(6) 

) 

IHR 

m 

JINT ( 

INDATA 

(7) 

) 

MIN 

m 

JINT ( 

INDATA 

(8) 

) 

11 

m 

JINT ( 

INDATA 

(9) 

) 

F10 

m 

INDATA 

(10) 
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F10B - INDATA (11) 
61 - INDATA (12) 


CALL TME ( MN . IDA . lYR , IHR . MIN . XUT . XLNG . SDA , 

SHA . DD . DY ) 

CALL TINF ( F10 . F10B , 01 , XLAT . SDA . SHA . DY , II . TE ) 

CALL JAC ( Z . TE , T2 . A(1) . A(2) . A(3) . A(4) . A(5) , A(6) . 
EM . DENS . DL ) 


DENL6 - 0. 

DUMkIY - DL 
DEN - DL 

IF ( Z .LE. 170. ) THEN 

CALL SLV ( DUWMY , Z . XUT , DD ) 
DENL6 - DUyMY 

END IF 


C 

C** ‘Fair* helium number density between base fairing height ( BFH ) and 500 km 
C 


IF ( Z. GE. 500. ) THEN 

CALL SLVH ( DEN . A(5) , XUT . SDA ) 

DL - DEN 

ELSE IF ( Z .GT. BFH ) THEN 

DHEL1 - A ( 5 ) 

DHEL2 - A ( 5 ) 

DL61 - DL 
DL62 - DL 

CALL SLVH ( DL62 , DHEL2 . XUT . SDA ) 

IH - Z 

CALL FAIRS ( DHEL1 , DHEL2 . DL61 . DLG2 . IH . FDHEL , FDL6 ) 
DL - FDLG 
A ( 5 ) - FDHEL 

END IF 


DL - DL + DENLG 

DENS - ie.**DL 

XUT - XUT e 57.29577951 


C Fill OUTDATA array 

OUTDATA (1) - TE 
OUTDATA (2) - TZ 

DO 80 I - 1 , 6 

OUTDATA (1+2) - 1.E6 e ( 10. •• A(I) ) 
80 CONTINUE 

OUTDATA (9) - EM 
OUTDATA 00) * dens • 1900. 

OUTDATA (11) - DL 
P - OUTDATA (10) e RGAS • TZ / EM 
OUTDATA (12) - P 

RETURN 

END 
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SUBRCXiTINE TME ( MN , IDA » lYR , IHR , MIN , XLAT , XLNG , 
SDA . SHA . DD , DY ) 


C**»*******«**«***«*««««*«**«********»*#****«****«***o****«***«*«>«**«***«****« 

C*» Subroutine *TME* performs the colculations of the solar declination * 

ongle and solar hour angle, • 

0*9 • 

C«* INPUTS: MN » month • 

C** IDA ■ doy • 

C** lYR - year • 

C*^ IHR - hour ♦ 

C** MIN » minute • 

C** XMJD«> mean Julian date * 

C** XLAT- latitude ( input-geocentric latitude ) • 

C** XLNG- longitude { input-geocentric longitude, -180, +180 ) • 

C*« • 

C** OUTPUTS: SDA - solar declination angle (rad) • 

C** SHA - solar hour angle (rod) • 

DD — day number from 1 JAN. • 

O** DY - DD / tropical year • 

C** Modified by Mike Hickey, USRA • 

DIMENSION IDAY(12) 

DATA IDAY / 31,28,31 ,30,31,30 ,31,31,30 ,31,30,31 / 

PARAMETER YEAR - 365.2422 

PARAMETER A1 - 99.6909833 , A2 - 36000.76892 

PARAMETER A3 - 0.00038708 , A4 - 0.250684477 

PARAMETER B1 - 0.0172028 , B2 - 0.0335 , B3 - 1.407 

PARAMETER PI - 3.14159265 , TPI - 6.28318531 

PARAMETER PI2 - 1.57079633 , PI32 - 4.71238898 

PARAMETER RAD^DEG - 0.017453293 

XLAT - XLAT / 57.29577951 
YR - lYR 

IF ( JM0D(IYR,4) .EQ. 0 ) THEN 

IF ( JMOD(IYR,100) .NE. 0 ) IDAY(2) - 29 I Century not a leap year 
ELSE 

IDAY(2) - 28 
END IF 

ID - 0 

IF ( MN. GT. 1 ) THEN 
DO 20 I - 1 , MH-1 
ID - ID + IDAY(I) 

20 CONTINUE 

END IF 

ID - ID + IDA 
DD - ID 
DY - DDAEAR 
C 

C«* Compute mean Julian date 
C 

XMJD - 2415020. + 365. • ( YR - 1900. ) + DD 

+ FLOATJ ( ( lYR - 1901 ) / 4 ) 

C 

C*« Compute Greenwich mean time in minutes GMT 
C 

XHR - IHR 
XMIN - MIN 

GMT - 60 e XHR + XMIN 

FMJD - XMJD - 2435839. + GMT / 1440. 

C 

C** Compute Greenwich mean position - GP ( in rad ) 

C 

XJ - ( XMJD - 2415020.5 ) / ( 36525.0 ) 

GP - AMOD ( A1 + A2 • XJ + A3 • XJ • XJ + A4 • GMT , 360. ) 

C 

C** Compute right ascension point - RAP ( in rad ) 

C 

C»* 1st convert geocentric longitude to deg longitude — west neg , + east 


IF ( XLNO ,CT, 180. ) XING - XLNG - 360. 
RAP - AMOD ( GP + XLNG , 360. ) 


C 

Compute celestial longitude - XLS ( in rad ) ~ zero to 2PI 

C 

Y1 - B1 • FMJD 

Y2 - 0.017202 • ( FMJD - 3. ) 

XLS - AMOD ( Y1 + B2 • SIN(Y2) - B3 , TPI ) 

C 

Compute solar declination angle SDA ( In rad ) 

C 

B4 - RAD.DEG • ( 23.4523 - 0.013 • XJ ) 

SDA - ASIN ( SIN ( XLS ) • SIN ( B4 ) ) 

C 

C** Compute right ascension of Sun - RAS ( in rad ) - - zero to 2PI 
C 

RAS - ASIN { TAN ( SDA ) / TAN ( B4 ) ) 

C 

C*« Put RAS in same quadrant os XLS 
C 

RAS - ABS ( RAS ) 

TEMP - ABS ( XLS ) 

IF ( TEMP. LE. PI .AND. TEMP.GT.PI2 ) THEN 
RAS * PI ~ RAS 

ELSE IF ( TEMP.LE.PI32 .AND. TEMP. GT. PI ) THEN 

ELSE IF OEMP^GT^Pm ) THEN 
RAS - TPI - RAS 

END IF 

IF ( XLS. LT. 0. ) RAS - -RAS 
C 

C** Compute solar hour angle — SHA ( in deg ) — * 

C 

SHA - RAP • RAD^DEG - RAS 

IF ( SHA.CT.PI ) SHA - SHA - TPI 

IF ( SHA.LT.-PI ) SHA - SHA + TPI 


RETURN 

END 
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SUBROUTINE TINF ( F10 . F10B . GI. XLAT, SDA . SHA . OY , 11 . TE ) 




C** Subroutine *TINF* calculates the exospheric temperoture according to 
C** L. Jocchia SAO 313, 1970 

C** 

C*e F10 — solar radio noise flux ( x E-22 Watts / m2 ) 

C*e F10B- 162-day averoge F10 

61 " geomagnetic activity index 

lAT « geographic latitude at perigee ( in rad ) 

C** SDA • solar declination ongle ( in rad ) 

C*e SHA « solar hour angle 

Ce* DY - D / Y ( day number / tropical year ) : 1 

Gee II ■ geomagnetic equation index ( 1 — GI»KP , 2— GI-AP ) 

Gee RE - diurnal foctor KP, F10B, AVG 

Ge* 

G*e CONSTANTS — C * solar activity variation 

Ce* — BETA , etc - diurnol variation 

G«e — D >■ geomagnetic variation 

C*« — E semiannual variation 

Ce* 

C«* Modified by Mike Hickey, USRA 

C«e*«**««*««e«*««eeeee***«*e«*«*****ee*e*eee«**ee«e«*e*e**«*«eee*e*e*ee«***»< 

PARAMETER PI - 3.14159265 , TPI - 6.28318531 
PARAMETER XM - 2.5 . XNN - 3.0 
C 

C*« Ci are solar activity variation vorlobles 
C 

PARAMETER Cl - 383.0 , C2 - 3.32 , C3 « 1.80 
C 

C** Di are geomagnetic variation variables 
C 

PARAMETER DI - 28.0 , D2 - 0.03 , D3 - 1.0 , D4 - 100.0 , D5 - -0.08 
C 

C«e Ei are semiannual variation variables 
C 

PARAMETER El - 2.41 , E2 - 0.349 , E3 - 0.206 , E4 - 6.2831853 
PARAMETER E5 - 3.9531708 . E6 - 12.5663706 , E7 - 4.3214352 
PARAMETER E8 - 0.1145 , E9 - 0.5 , E10 - 6.2831853 
PARAMETER Ell « 5.9742620 , E12 - 2.16 

PARAMETER BETA - -0.6457718 , GAMMA - 0.7504916 , P - 0.1047198 
PARAMETER RE - 0.31 
C 

Ge* solar activity variation 
C 

TC - Cl + C2 • F10B + C3 ♦ ( F10 - F10B ) 

C 

C*e diurnal variation 
C 

ETA - 0.5 • ABS ( XLAT - SDA ) 

THETA - 0.5 • ABS ( XLAT + SDA ) 

TAU » SHA 4 BETA 4 P • SIN ( SHA 4 GAMMA ) 

IF ( TAU. GT. PI ) TAU - TAU - TPI 
IF ( TAU. LT.-PI ) TAU - TAU 4 TPI 

A1 - ( SIN ( THETA ) )**XM 

A2 - ( COS ( ETA ) )**XM 

A3 - ( COS ( TAU / 2. ) )**XNN 

B1 - 1.0 4 RE e A1 

B2 - ( A2 - A1 ) / B1 

TV - B1 • ( 1 . 4 RE e B2 * A3 ) 

TL - TC • TV 
C 

Ce* geomagnetic variation 
C 

IF ( I1.EQ.1 ) THEN 

TG - DI • GI 4 D2 • EXP(GI) 

ELSE 

TG - D3 • GI 4 D4 ♦ ( 1 - EXP ( D5 • GI ) ) 

END IF 
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C*« semiannual variation 
C 

G3 - 0.5 • ( 1.0 + SIN ( E10 • DY + Ell ) ) 

63 - 63 E12 

TAU1 - DY + E8 • ( 63 - E9 ) 

61 - E2 + E3 • ( SIN ( E4 • TAU1 + E5 ) ) 

62 - SIN ( E6 • TAU1+ E7 ) 

TS - El + F10B • 61 • 62 

C 

exospheric temperature 
C 

TE - TL + T6 + TS 

RETURN 

END 


SUBROUTINE JAC ( Z , T , TZ . AN , A02 , AO . AA . AHE . AH , EM , 
DENS , DL ) 




c** 

C»* Subroutine *JAC’ calculates the temperature TZ , the total density DENS 
C** and its logarithm DL, the mean molecular weight EM, the individual 
C** specie number densities for N, 02, 0, A, HE and H ( each preceded with 
C*» an ’A’ ) at altitude Z given the exospheric temperature T. 

C»* This subroutine uses the subroutine 'GAUSS* and the function 
C*e subprogroms 'TEMP' and *MOU_WT'. 

C*e 

C*e Rewritten by Mike Hickey, USRA 

* 


DIMENSION ALPHA(6) , EI(6) , DI(6) . DIT(6) 

REAL*4 MOL_WT 

PARAMETER AV - 6.e2257E23 
PARAMETER ON - .78116 
PARAMETER 002 - .26955 
PARAMETER OA - .669343 
PARAMETER OHE - 1.289E-65 
PARAMETER RGAS - 8.31432 
PARAMETER PI - 3.14159265 
PARAMETER T6 - 183. 

GRAVITY ( ALTITUDE ) - 9.86665 / ( ( 1. + ALTITUDE / 6.356766E3 )**2 ) 


DATA ALPHA / 6.6 , 6.6 , 6.6 , 6.6 , -.386 , 6.6 / 

DATA El / 28.6134 , 31.9988 . 15.9994 , 39.948 , 4.6626 , 1.66797 / 


TX - 444.3867 + .62385 e T - 392.8292 e EXP ( -.6621357 e T ) 

A2 - 2. e (T-TX) / PI 

TX_T6 - TX - T6 

T1 - 1 .9 e TX_T6 / 35. 

T3 - -1.7 e TX_T6 / ( 35.e.3 ) 

T4 - -6.8 e TX_T6 / ( 35. **4 ) 

TZ - TEMP ( Z . TX , T1 . T3 . T4 . A2 ) 


C*e SECTION 1 

C*e 

A - 96. 

D - AMIN1 ( Z , 16S. ) 

C Integrate gM/T from 96 to minimum of Z or 165 km :- 
CALL GAUSS ( A, D. 1 , R. TX , T1 , T3 , T4 , A2 ) 

C The number 2.1926E— 8 > density x temperature/mean molecular weight at 96 km. 
EM - MOU.WT ( D ) 

TD - TEMP ( D . TX , T1 . T3 , T4 , A2 ) 

DENS - 2.1926E-8 e eM e EXP( -R / RGAS ) / TD 

FACTOR - AV e DENS 
PAR - FACTOR / EM 
FACTOR - FACTOR / 28.96 


C For altitudes below and at 165 km calculate the individual specie number 
C densities from the mean molecular weight and total density. 

IF ( Z. LE. 165 ) THEN 

DL - AL0G16 ( DENS ) 

AN - AL0G16 ( ON e FACTOR ) 
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AA - ALOG10 ( QA • FACTOR ) 

AHE - ALOG10 ( QHE • FACTOR ) 

AO - ALOG10 ( 2. • PAR # ( 1 --EM / 28.96 ) ) 

A02 - ALOG10 ( PAR • ( EM • ( 1.4002 ) / 28.96-1. ) ) 
AH - 0. 

C 

C«* Return to calling program 
C 

RETURN 


END IF 


C*« SECTION 2 : This section is only performed for altitudes above 105 km 

C** 


Note that having reached this section means that D in section 1 is 105 km. 


C Calculate individual specie number densities from the total density and mean 
C molecular weight at 105 km altitude. 

DI(1) - ON • FACTOR 

DI(2) - PAR e (EM • (1.4002) / 28.96-1.) 

DI(3) - 2. • PAR • (1.- EM / 28.96) 

DI(4) - QA • FACTOR 
DI(5) - QHE • FACTOR 

C Integrate g/T from 105 km to Z km :- 

CALL GAUSS ( D, Z. 2, R. TX , T1 , T3 , T4 . A2 ) 

DO 41 1-1,5 

DIT(I) - DI(I) e ( TO / TZ ) •♦(1 .4-ALPHA(I)) • EXP( -EI(I) e R / RGAS) 
IF ( DIT(I), LE. 0. ) DIT(I) - 1.E-6 
41 CONTINUE 


C*« This section calculates atomic hydrogen densities above 500 km altitude. 
C»* Below this altitude , H densities are set to 10**-6. 

C*^ SECTION 3 

C** 

IF ( Z .GT. 500. ) THEN 
A1 - 500. 

S - TEMP ( A1 . TX . T1 , T3 . T4 . A2 ) 

DI(6) - 10. •• ( 73.13 - 39.4 • ALOG10 (S) 4- 5.5 • ALOG10(S) *ALOG10(S)) 
CALL GAUSS ( A1 , Z, 7, R. TX , T1 . T3 , T4 , A2 ) 

DIT(6) - 01(6) • (SA2) • DCP ( -EI(6) • R / RGAS ) 

ELSE 

DIT (6) - 1.0 

END IF 


C For altitudes greater than 105 km , calculate total density and mean 
C molecular weight from individual specie number densities. 

DENS-0 

DO 42 1-1,6 

DENS - DENS + EI(I) • DIT(I) / AV 

CONTINUE 
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EM - DENS • AV / ( DIT(1)+0IT(2)+DIT(3)4DIT(4)+01T(5)+0IT(6) ) 
OL - ALOGie (DENS) 

AN - AL06ie(DlT(1)) 

A02 - ALOGie(DIT(2)) 

AO - AL06ie(DIT(3)) 

AA - AL0G10(DIT(4)) 

AHE - ALOG10(DIT(5)) 

AH « AL06ie(DIT(6)) 


RETURN 



FUNCTION TEMP ( ALT , TX , T1 , T3 . T4 , A2 ) 


* 

c** 

C** Function subprogram *TEMP* calculates the temperature at oltitude ALT 
C** using equation (10) for altitudes between 90 and 125 km and equation 
C** (13) for altitudes greater than 125 km , from SAO Report 313. 

C** 

C** Written by Mike Hickey, USRA 

C***********«****************«********««**«***«*«**********«****«******«** 


parameter BB - 4.5E-6 

U - ALT - 125. 

IF ( U .GT . 0. ) THEN 

TEMP ■« TX + A2 • ATAN (T1eUe(1.+BBe (U**2.5)) / A2 ) 

ELSE 

TEI^ - TX + T1 • U + T3 • (U**3) + T4 • (U**4) 

END IF 
END 


REAL FUNCTION M0U_WT*4 ( A ) 


** 

Subroutine 'MOU.WT* calculates the molecular weight for attitudes •• 

between 90 and 105 km according to equation (1) of SAO report 313. •• 

C** Otherwise, MOU_WT is set to unity. •• 

•• 

Written by Mike Hickey, USRA •• 

C****************************************************************************** 


DIMENSION B (7) 

DATA B / 28.15204 , -0.085586, 1.284E-4. -1.0056E-5, -1.021E-5, 
1.5044E-6, 9.9826E-8 / 

IF ( A. GT. 105. ) THEN 

MOU_WT - 1 . 

ELSE 

U - A - 100. 

MOU-WT - B (1) 

DO 1 1-2,7 

MOU-WT - MOU.WT + B (I) • U ♦* ( 1-1 ) 

1 CONTINUE 

END IF 
END 
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SUBROUTINE GAUSS ( Z1 . Z2 . NMIN . R . TX . T1 , T3 . T4 . A2 ) 




C** Subdivide total integration— oltitude ronge into intervals suitoble for • 
C** applying Gaussian Quadrature , set the number of points for integration * 
C** for each sub— interval , and then perform Gaussian Quodrature. • 

C*e Written by Mike Hickey, USRA, NASA/MSFC. E044. July 1988. • 



REAL*4 ALTMIN (9) . C(8.6). X(8,6). MOL_WT 
INTEGER NG (8) , NGAUSS . NMIN . J 

GRAVITY ( ALTITUDE ) - 9.80665 / ( ( 1. + ALTITUDE / 6.356766E3 )**2 ) 

DATA ALTMIN / 90.. 105., 125., 160., 200., 300., 500., 1500., 2500. / 
DATA NG / 4. 5, 6. 6. 6. 6, 6. 6/ 

C Coefficients for Gaussian Quadrature ... 


DATA C / 

.5555556 , 

.8888889 , 

.5555556 , 

0000000 . 1 

n«3 



.0000000 , 

.0000000 . 

.0000000 , 

0000000 , 1 

n«3 



.3478548 • 

.6521452 . 

.6521452 . 

3478548 , 1 

n»4 



.0000000 , 

.0000000 , 

.0000000 . 

0000000 . 1 

n-4 



.2369269 , 

.4786287 . 

.5688889 . 

4786287 , I 

n*5 



.2369269 . 

.0000000 , 

.0000000 . 

0000000 , 1 

n«^ 



.1713245 . 

.3607616 . 

.4679139 , 

4679139 , ! 

n*“6 



.3607616 , 

.1713245 . 

.0000000 , 

0000000 • 1 

n— 6 



.1294850 , 

.2797054 * 

.3818301 . 

4179592 , 1 

n«7 



.3818301 , 

.2797054 , 

.1294850 , 

.0000000 , 1 

n-7 



.1012285 , 

.2223810 , 

.3137067 , 

3626838 , ! 

n-8 


.3626838 . .3137067 , 

C Abscissas for Gaussian Quadrature ... 

.2223810 . 

1012285 / 1 

n»8 


DATA X / 

-.7745967 

, .0000000 

, .7745967 

. .0000000 

. 1 

n-3 


.0000000 

, .0000000 

, .0000000 

, .0000000 

. 1 

n*3 


-.8611363 

. -.3399810 

, .3399810 

, .8611363 

• 1 

n— 4 


.0000000 

» .0000000 

, .0000000 

, .0000000 

0 1 

n— 4 


-.9061798 

. -.5384693 

, .0000000 

, .5384693 

. ! 

n-5 


.9061798 

, .0000000 

. .0000000 

, .0000000 

. 1 

n*5 


-.9324695 

, -.6612094 

. -.2386192 

e .2386192 

• I 

n»6 


.6612094 

. .9324695 

, .0000000 

, .0000000 

• 1 

n*6 


-.9491079 

, -.7415312 

. -.4058452 

, .0000000 

• ! 

n-7 


.4058452 

, .7415312 

, .9491079 

, .0000000 

, 1 

n— 7 


-.9602899 

, -.7966665 

, -.5255324 

, -.1834346 

. 1 

n-8 


. 1834346 

, .5255324 

, .7966665 

, .9602899 / 1 

n— 8 


R - 0.0 


DO 2 K - NMIN . 8 

NGAUSS - NG (K) 

A - ALTMIN (K) 

D - AMIN1 ( Z2 . ALTMIN (K-f1) ) 

RR - 0.0 

DEL - 0.5 • ( D - A ) 

J - NGAUSS - 2 

DO 1 I - 1 , NGAUSS 

Z - DEL e ( X(I,J) + 1. ) + A 

RR - RR + C(I.J) • MOU_WT(Z) e GRAVITY(Z) / TEMP ( Z.TX.T1 .T3.T4.A2 ) 


1 CONTINUE 

RR - DEL • RR 
R - R + RR 

IF ( D .EQ. Z2 ) RETURN 

2 CONTINUE 



ORIGINAL PA.GE IS 
OF POOR QUALITY 
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RETURN 

END 



SUBROUTINE SLV ( DEN , ALT , XLAT , DAY ) 


C* 

C* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


Subroutine *SLV* computes the seasonal— I at i tudi nol variation of density 
in the lower thermosphere in accordance with L« Jacchia» SAO 332» 1971. 
This affects the densities between 90 and 170 km. This subroutine need 
not be called for densities above 170 km, because no effect is observed. 

The variation should be computed after the calculation of density due to 
temperature variations and the density ( DEN ) must be in the form of o 
base 10 1 09 . No adjustments are made to the temperature or constituent 
number densities In the region affected by this variation. 

DEN « density (Iog 10 ) 

ALT - altitude (km) 

XUT - lotitude (rad) 

DAY m day number 


initialize density (DEN) « 0.0 
C 

DEN - 0.0 
C 


check if altitude exceeds 170 km 
C 

IF ( ALT. GT. 170. ) RETURN 
C 


C** compute density change In lower thermosphere 
C 

Z - ALT - 90. 

X - -0.0013 • Z • Z 
Y - 0.0172 • DAY + 1.72 
P - SIN (Y) 

SP - ( SIN (XLAT) ) •♦2 
S - 0.014 • Z • EXP (X) 

D - S e P e SP 
C 

check to compute absolute value of *XLAT* 

C 

IF ( XLAT. LT. 0. ) D - -0 
DEN - D 


RETURN 

END 


SUBROUTINE SLVH ( DEN . OENHE . XLAT , SOA ) 


C** Subroutine ‘SLYH* computes the seasonal-latitudinal varaltion 
C** helium number density according to L. Jacchia, SAO 332, 1971. 
correction is not important below obout 500 km. 






DEN - 

density (Iog10) 

c** 

DENHE - 

helium number density (Iog10) 


XLAT - 

latitude (rod) 


SDA - 

solar declination angle (rod) 


D0 - 10. DENHE 



A - ABS ( 0.65 • ( SDA / 0.40909079 ) ) 

B - 0.5 e XLAT 
C 

C** Check to compute absolute value of *B* 

C 

IF ( SDA. LT. 0. ) B - -B 
C 

C*« compute X, Y. DHE ond DENHE 
C 

X « 0.7854 - B 
Y - ( SIN (X) ) •• 3 
DHE- A • ( Y - 0.35356 ) 

DENHE - DENHE + DHE 
C 

C** compute helium number density change 
C 

D1 - 10. *« DENHE 

DEL- D1 - D0 

RHO- 10. ** DEN 

DRHO - ( 6.646E-24 ) • DEL 

RHO - RHO 4- DRHO 

DEN - ALOG10 (RHO) 

RETURN 

END 


of the 
This 
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SUBROUTINE FAIRS ( DHEL1 .DHEL2 ,DLC1 ,DLG2 ,IH ,FDHEL ,FDLG ) 


C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

Ct 

C* 

C* 

C* 

C* 


Thi« subroutine fairs between the region above 500 km, which invokes the 
seosona l-iat i tudi nal variation of the helium number density ( subroutine 
SLVH ), and the region below, which does not invoke any seasonal- 
latitudinal variation at all. 

INPUTS: DHEL1 - helium number density before invoking SLVH 

DHEL2 - helium number density after invoking SLVH 
DLG1 - total density before invoking SLVH 
DLG2 « total density after invoking SLVH 
IH - height ( km ) — INTEGER 
IBFH - base fa i ring height ( km ) — INTEGER 
OUTPUTS: FDHEL • foired helium number density 
FDLG « faired total density 

Written by Bill Jeffries, CSC. Huntsville, AL. 
ph. (205) 830-1000, x311 


DIMENSION CZ ( 6 ) 

DATA CZ / 1.0, 0.9045085, 0.6545085, 0.3454915, 0.0954915, 0.0 / 
PARAMETER IBFH - 440 

C Height index 

I - ( IH - IBFH ) /10 + 1 
C Non-SLVH fairing coefficient 
CZI - CZ ( I ) 

C SLVH fairing coefficient 
SZI - 1.0 - CZI 
C Fai red dens! ty 

FDLG - ( DLG1 e CZI ) + ( DLG2 • SZI ) 

C Foired helium number density 

FDHEL - ( DHEL1 • CZI ) + ( DHEU e SZI ) 

RETURN 

END 


☆U.S. GOVERNMENT PRINTING OfTICE: WM - 530-050/10107 
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